module testmoments2
implicit none
contains
function emempkememp(b,eta,si)
!E(kron(emem',emem'))
use normalmoms
use matrixfuns
use linear_operators
double precision eta,si
double precision b(:)&
&,emempkememp(1:size(b)*size(b),1:size(b)*size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))&
&,vb(1:size(b),1:1)
double precision h1,h2,h3,h4,p1,p2

call paramstest2(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
p1=h4-4*h3+6*h2-4*h1+1 !E((h-1)^4)
p2=h3-2*h2+h1 !E(h*(h-1)^2)
!!
vb(:,1)=b
emempkememp=c**4.0*p1*kron(vb.xt.vb,vb.xt.vb)&
          &+c**2.0*p2*kron(vb.xt.vb,am.xt.am)&
		  &+c**2.0*p2*(vec(vb.xt.vb).xt.vec(am.xt.am))&
		  &+c**2.0*p2*(kron(vb,am).xt.kron(am,vb))&
		  &+c**2.0*p2*(kron(am,vb).xt.kron(vb,am))&
		  &+c**2.0*p2*(vec(am.xt.am).xt.vec(vb.xt.vb))&
		  &+c**2.0*p2*kron(am.xt.am,vb.xt.vb)&
		  &+h2*(((eye(size(b)**2)+commatrix(size(b),size(b))).x.kron(am.xt.am,am.xt.am))&
		      &+(vec(am.xt.am).xt.vec(am.xt.am)))
end function emempkememp
!!
function emkememp(b,eta,si)
!E(kron(em,em)*em')
use normalmoms
use matrixfuns
use linear_operators
double precision eta,si
double precision b(:)&
&,emkememp(1:size(b)*size(b),1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,p1,p2

nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si

h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
p1=h3-3.0*h2+3.0*h1-1 !E((h-1)^3)
p2=h2-h1 !E(h*(h-1))



call paramstest2(b,am,c,nu,gam,eta,si)
emkememp=emkeheips(b,am)

	

contains 
function emkeheips(b,am)
!E(kron(em,eh)*ei')
use normalmoms
use matrixfuns
use linear_operators
double precision b(:),emkeheips(1:size(b)*size(b),1:size(b))
double precision am(1:size(b),1:size(b))&
&,vb(1:size(b),1:1)

vb(:,1)=b

emkeheips=c**3.0*p1*(kron(vb,vb).xt.vb)&
        &+c*p2*((kron(vb,am)+kron(am,vb)).xt.am)&
		&+c*p2*(vec(am.xt.am).xt.vb)
end function emkeheips
end function emkememp
!!
function jimemkememp(b,eta,si)
!E(ji_m*kron(em,em)*em')
use normalmoms
use matrixfuns
use linear_operators
double precision eta,si
double precision b(:),jimemkememp(1:size(b)*size(b),1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,p1,p2,p3
integer ind
nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si

h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
p1=h5-5*h4+10*h3-10*h2+5*h1-1  !E((h-1)^5)
p2=h4-3*h3+3*h2-h1  !E(h(h-1)^3)
p3=h3-h2  !E(h^2(h-1))

call paramstest2(b,am,c,nu,gam,eta,si)
jimemkememp=jijemkeheips(b,am)
contains 
function jijemkeheips(b,am)
!E(ji_jkron(em,eh)*ei')
use normalmoms
use matrixfuns
use linear_operators
double precision b(:),jijemkeheips(1:size(b)*size(b),1:size(b))
double precision am(1:size(b),1:size(b))&
&,vb(1:size(b),1:1)
integer nk

vb(:,1)=b

nk=size(b)


jijemkeheips=c**5.0*p1*dot_product(b,b)*(kron(vb,vb).xt.vb)&
           &+c**3.0*p2*quadnorm1(am.tx.am)*(kron(vb,vb).xt.vb)&
		   &+2.0*c**3.0*p2*((kron(vb,vb).xt.vb).x.(am.xt.am))&
		   &+2.0*c**3.0*p2*(((kron(vb,am)+kron(am,vb)).xt.am).x.(vb.xt.vb))&
		   &+c**3.0*p2*dot_product(b,b)*((kron(vb,am)+kron(am,vb)).xt.am)&
		   &+c*p3*((kron(vb,am)+kron(am,vb)).x.(quadvec1(am.tx.am).xt.am))&
		   &+c**3.0*p2*dot_product(b,b)*(vec(am.xt.am).xt.vb)&
		   &+c*p3*((kron(am,am).x.vec(quadvec1(am.tx.am))).xt.vb)&
		   &+2.0*c*p3*(kron(am,am)&
		   &.x.((eye(nk**2)+commatrix(nk,nk)+(vec(dble(eye(nk))).xt.vec(dble(eye(nk)))))&
		   &.x.(kron(am.tx.vb,dble(eye(nk))).xt.am)))

end function jijemkeheips
end function jimemkememp


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine paramstest2(b,am,c,nu,gam,eta,si)
use linear_operators
use bessels
use ghs
double precision, intent(in):: b(:),eta,si
double precision, intent(out):: nu,gam,c,am(1:size(b),1:size(b))
double precision a,normb
integer k,n

nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si
c=ghcfun0(b,eta,si)
if (dot_product(b,b).gt.1.0d-3) then
	am=eye(size(b))+((c-1)/dot_product(b,b))*(matvec(b) .xt.matvec(b))
	am=.t.chol(am)
else 
	a=besseld(nu+1.0d0,gam)-1.0d0
	normb=dot_product(b,b)
	am=eye(size(b))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(b) .xt. matvec(b))
	am=.t.chol(am)
end if

end subroutine paramstest2
function eh(nu,gam,k)
!E(ht**k)
use bessels
double precision nu,gam,eh
integer k,i
eh=1
do i=1,k
	eh=eh*besselr(nu+i-1,gam)/besselr(nu,gam)
end do

end function eh
end module testmoments2